
function [p q] = SemiImplicit_BDF(alpha,beta)
global U  h  g  p  q  dt nt dx x h1 p_global;

dx1=1/dx;
dx2 = 1/(dx*dx);
dt2= 1/dt;

%---- Construct the spatial Discretization Matrix------%
A = zeros(x,x);
rhs = zeros(x,1);

A(1,x)=-U*0.5*dx1;
A(1,2) = U*0.5*dx1;
for k=2:x-1
    A(k,k-1)=-U*0.5*dx1;
    A(k,k+1) = U*0.5*dx1;
    %    A(k,k)=1/dt;
end
A(x,x-1)=-U*0.5*dx1;
A(x,1)=U*0.5*dx1;


B=zeros(x,x);

B(1,x)=h*dx2;
B(1,2)=h*dx2;
B(1,1)=-2*h*dx2;
for k=2:x-1
    B(k,k)= -2*h*dx2;
    B(k,k-1)=h*dx2;
    B(k,k+1)=h*dx2;
end
B(x,1)=h*dx2;
B(x,x-1)=h*dx2;
B(x,x)=-2*h*dx2;

C = g*eye(x);
I=eye(x);
itermax=10;
beta1=0.5;
D = (dt2*I + beta1*A);
E = (dt2*I + beta*A);

% Initiate Time Loop
for n=2:nt+1;
    % For the first time steps, implement the First Order methods
    if (n==2)
        for iter =1:itermax+5
            % advance q explicitly by Forward Euler method
            if (iter==1)
                qex =q;
                %qex(1:x,1) = -dt*(C*p(1:x,1)+A*q(1:x,1)) + q(1:x,1);
            end
            % solve p implicitly Trapezoidal
            rhs = dt2*p(1:x) - (1-beta1)*(A*p(1:x,1)+B*q(1:x,1)) - beta1*B*qex(1:x,1);
            p_np1 =D\rhs;
            % Correct q implicitly Trapezoidal
            rhs = dt2*q(1:x) - (1-beta1)*(C*p(1:x,1)+A*q(1:x,1)) - beta1*C*p_np1(1:x,1);
            qex  =  D\rhs;
        end
        q_nm1 = q;
        p_nm1 = p;
        p_n   = p_np1;
        q_n   = qex;
        
    else
        % For the second time step onwards, second order methods can be
        % implemented
        for iter =1:itermax
            % advance q explicitly by Leapfrog
            if(iter==1)
                qex =q_n;
                %qex(1:x,1) = -2*dt*(C*p_n(1:x,1)+A*q_n(1:x,1)) + q_nm1(1:x,1);
            end
            
            % Solve for p implicitly
            rhs = (dt2*(1+alpha))*p_n(1:x,1) +(alpha-(1-beta))*(A*p_n(1:x,1)+ B*q_n(1:x,1)) -beta*(B*qex(1:x,1)) -(alpha*dt2)*p_nm1(1:x,1);
            p_np1=E\rhs;
            % Correct q implicitly
            
            rhs = (dt2*(1+alpha))*q_n(1:x) + (alpha-(1-beta))*(C*p_n(1:x,1)+A*q_n(1:x,1)) -beta*(C*p_np1(1:x,1)) -(alpha*dt2)*q_nm1(1:x,1);
            qex=E\rhs;
        end
        q_nm1 = q_n;
        p_nm1 = p_n;
        p_n   = p_np1;
        q_n   = qex;
        p = p_n;
    end
    time =n*dt;
    if rem(time,5)==0
        k=time/5;
        p_global(:,k) = p;
        refreshdata(h1,'caller') % Evaluate p in the function workspace
        drawnow
    end
end
display('Completed Successfully');